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ABSTRACT 

Observations show that there is a proton spectral “break” with B[,reak at 1-lOMeV 
in some large CME-driven shocks. Theoretical model usually attribute this phe¬ 
nomenon to a diffusive shock acceleration. However, the underlying physics of the 
shock acceleration still remains uncertain. Although previous numerical models can 
hardly predict this “break” due to either high computational expense or shortcomings 
of current models, the present paper focuses on simulating this energy spectrum in 
converged two shocks by Monte Carlo numerical method. Considering the Dec 13 
2006 CME-driven shock interaction with an Earth bow shock, we examine whether 
the energy spectral “break” could occur on an interaction between two shocks. As re¬ 
sult, we indeed obtain the maximum proton energy up to lOMeV, which is the premise 
to investigate the existence of the energy spectral “break”. Unexpectedly, we further 
find a proton spectral “break” appears distinctly at the energy ~5MeV. 

Subject headings: acceleration of particles-methods:numerical-shock waves-solar 
wind-Sun:coronal mass ejections(CMEs) 
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1. Introduction 


Diffusive shock acceleration (DSA) is regarded as the most efficient mechanism for galactic 
cosmic rays (CR). In the recent years several in situ observations indicate the supernova rem¬ 
nants (SNRs) are the essential candidates for the sources of the CR’s proton spectrum up to the 
“knee” at a few PeV. Since the TeV y-ray from Crab Nebula were first clearly detected by imag¬ 
ing air Cherenkov telescope (lACT) in 1989, lACTs have been extensively constructed and are 
operating around the wo rld. There are about one hundred y-ray sources are identified up to now 
(lAmenomori et alJl2009t> . Among of those sources, IC443 and W44 are the highest-significance 
SNRs in the second Fermi-LAT suited for a detailed study of their y-ray spectra. IC443 and W44 
are located at distances of 1.5kpc and 2.9kpc, respectively. Their GeV y-ray spectra with high- 
energy breaks at 20GeV and 2GeV for IC443 and W44 can trace the parent protons escaped from 
SNRs shock and collided with the nearby molecular clouds (MCs). In addition, both the y-ray 
emissions with a low-energy break at ~ 200MeV also can be inte rpreted by effect of pion bump 
between SNRs protons and the MCs (i.e. p-i-p^He-i-7t®, 7i®^2y) (lAckermann et al.ll2013b . This 
y-ray luminosity model can be used to explain the CRs energy spectral “break”. 

In the interplanetary (IP) space, there are also find the similar energy spectral “break” in the 
IP shocks. There are six events with hard energy spectra occurred on Nov. 6, 97, Feb. 15, 01, Jan. 
20, 05, Sep. 7, 05, Dec. 5, 06, an d Dec. 13,06. These six large events all have spectral “breaks” at 
the energy range of ~l-10MeV (iMewaldt et al.ll2008ll . In addition, there are the six largest events 
on the solar cycle 23 list as Jul. 14, 00, Nov. 8, 00, Sep. 24, 01, Nov. 04, 01, Nov. 22, 01, Oct. 
28, 03. These events all have spectra that roll over in similar fashion beyond ~50MeV. In present 
paper, we discuss the Dec. 13, 06 event, which proton fluxes are based on spacecraft of ACE, 
STEREO, and SAMPEX. 

Although a number of in situ observations exhibit the CR’s proton spectral “break” associated 
with either galaxy source or solar source, there is still no reliable prediction of this “break” by 
numerical methods. According to the theoretic model of DSA, CR’s power-law spectrum span a 
very large energy range from IKeV to a few 100EeV(~ lO^^eV). If one plans to simulate the total 
CR’s energy spectrum, it’s very hard for performing so expensive computer program. In addition, 
numerical simulation is usual to built a simple DSA model with a short size of the diffusive region 
ahead of shock. If the energy spectral “break” associated with a large diffusive size, the simple 
numerical models would hardly include this energy spectral “b reak” in their sirnulation results. 
Monte Carlo (MC) method can easily treat thermal ion injection (IWang et al.ll201lL 120131) . In MC 
method, the scattering mean free path is assumed to be some function of the particle rigidity, so 
this treatment is able to follow the evolution of individual ions long enough to model acceleration 
to high energies. However, the acceleration efficiency, as well as the maximum particle energy, are 
limited by the finite size of the accelerating region as parameterized by the free escape boundary 
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(FEB). I Vladimirov et alJ (l2008r) describes the escape of particles from the SNR shock with an 
assumed FEB far upstream of the shock. The distance to the FEB is a free parameter of the 
simulation that con trols the maximum energy of accelerated particles and the escaping CR flux. 
Ellison et al.l (119901) presented an ion spectra with a m aximum particle en ergy no more than IMeV 
by applying an fixed FEB ahead of the bow shock. iKnerr et al.l (119961) and IWang et al. I (120121) 
forward the maximu r n part icle energy achieving to ~4MeV by applying a moving FEB with the 
shock. I Wang et al. I (12015b investigated that the maximum particle energy limited in FEB can 
achieve a saturation at ~5.5MeV using different scattering mean free path functions. 

Since the cosmic rays are important both dynamically and diagnostically it is essential that we 
understand their acceleration, transport, radiative emissions, and interaction with other components 
of astrophysical environments. In particular, the CRs spectrum shape is usually referred to as a 
knee-ankle structure with the “knee” at a few PeV and the “ankle” at a few EeV. Similarly, the 
IP shock energy spectrum with a “break” at a few MeV. There are some debated understandings 
for the energy spectral “break”. One of these understandings has ever proposed that the “break” 
determined by the Larmor radius for ions. A heavy nucleus with charge Z has maximum energies 
Z times higher than a proton with same Lar mor radius . However, the CRs consist mainly of 
protons with smaller numbers of other nuclei (lBellll2013b . Another point argues that this “break” 
can probably be associated with the leakage mechanism . The CRs drive Alfven waves efficiently 
enough to build a tr ansport barrier that strongly reduces the leakage leading a spectral “break” 
(IMalkov et al .112013b . Although the turbulent field exist a limited boundary, it would not be the 
main reason for producing the energy spectral “break”. According to the DSA, the turbulent field in 
the precursor is generated by the gradient of CRs pressure, which is faded gradually but not steeply. 
However, if a SNR shock approaches a MCs or a pre-supemova swept-up shell with a significant 
amount of neutrals, confineme nt of accelerated pa rticles deteriorates. This would lead to a energy 
spectral “bre ak” (see review i n 
proposed by ISchneiden (11993b : 


Bvkov et al.ll2013b. In a ddition, there is a multiple shocks model 


Melrose & Pope.l (1 1993b . It is assumed that the medium is highly 


turbulent and that the number of shocks are propagating through it. CRs particles are injected 
into the system and accelerated by one or more shocks before they escaped from the system. This 
model may also used to study the particle spectral features such as “breaks”. In present paper, we 
propose that a new collided shocks model could probably inform the energy spectral “break” in the 
interacted precursor regions. 


How to directly follow the higher energy spectral “tail” for understanding the “break” is pos¬ 
ing a challenge to the numerical methods. Present work treats the interaction of CME-drive shock 
with the Earth bow shock to try investigate the energy spectral “break” as described in the ob¬ 
served Dec 13 2006 event. In converged two shocks, we apply the Monte Carlo simulation method 
without FEB. There are two reasons make it is possible to verify the spectral “break”. Firstly, the 
double shocks can provide enough energy injection to produce the required highest energy “tail”. 
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Secondly, due to no FEB, the accelerated particles can freely interact between double shocks and 
renew the standard DSA power-law shape. 

In present paper, we do simulations to further investigate the energy spectrum in converged 
two shocks by using Monte Carlo method. In section 2, we introduce the converged shocks model 
specifically. In section 3, the simulated results and analysis are presented. In the last section, we 
give summaries and conclusions. 


2. Model 


Halo CMEs were observed by the SOHO/LASCO coronagraphs in association with the events 
of 13 December 2006, with speeds of 1774kms“^. The flux spectra of protons in the SEP event by 
ACE, STEREO, and SAMPEX instruments show a “break” at ~l-10MeV. Spectra from GOES- 
11 also agree over the region from 5MeV to lOOMeV. Although the broken spectra would be little 
debated due to system errors from multiple spacecraft, it is hard to obtain a completely spectra for a 
large energy span from an individual spacecraft. So we hopefully do a simulation to obtain an entire 
spectra, which can cover the energy range from thermal energy below ~0.1MeV to superthermal 
energy beyond ~10MeV. 


Fig. [H shows a schematic diagram of the converged shocks model. The left reflective wall 
represents a CME, which produces a shock by a bulk flow speed of UO 2 . The right reflective wall 
represents the Earth, which informs a bow shock by an opposite bulk flow speed of mOi. Their 
relative speed between two bulk flows is equal to m = |m0i| -|- |m 02 |, which can equivalently be 
taken as the relative movement between two reflective walls with opposite velocities aOi, MO 2 in 
the laboratory reference frame. Similarly, we can take both downstream bulk flow speeds with 
the same velocity zero in the laboratory reference frame. This model describes the double shocks 
interaction occurred on the 13 December 2006 nearby Earth. Accordin g to the Wind m agnetic 
cloud list, the cloud axis direction is 6=27° ,(|)=85° in GSE coordinates (ILiu et al.ll2008b . In this 
model, we deflne the bulk flow direction to the interplanetary magnetic field (IMF) direction with 
an oblique factor of cos(0). So we take the relative speed between two bulk flows as value for 
~1600km“^ aligned to the IMF. Both two shocks are produced by the same bulk flow speed value 
for |m0i| = |M02|=~800km“\ but with opposite direction in the laboratory frame. Also both two 
reflective walls produce shocks propagating with opposite velocities of Vghi and v^/j 2 , respectively. 


In this Monte Carlo method, we apply an initial number density of particles no obeying a 
Maxwellian distribution with a random thermal velocity vq in the unshocked upstream region. 
Initial particles with their upstream bulk flow speeds move to their corresponding reflective walls 
in the both sides of the simulation box. Each bulk flow is reflected and forms the higher shocked 
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The Schematic Diagram of the Simulation Box in Double Shocks Model 



Fig. 1.— A schematic diagram of the simulation box. The left reflective wall represents a CME and 
produces one shock propagating from the left boundary to the middle of the simulation box. The right 
reflective wall represents the Earth and produces another shock involuting from the right boundary to the 
middle of the simulation box. Particles diffuse in the precursor regions between two shock fronts when they 
approach together near the middle of the simulation box. 

densities as downstream region in each side of the simulation box. With the downstream density 
achieving to a saturation, the shock fronts smoothly evolute forward to the middle of the box with 
shock velocity Vshi and Vshi, respectively. When the two shock fronts propagate more and more 
close, then the two precursors have an interaction with the time at the middle of the simulation box. 
It is this processes the particles are able to gain their additional energies by crossing two shock 
fronts. Unlike previous single shock model, FEB is not included in the front of the converged two 
shocks. Due to the two reflective walls in this model can play a role on preventing the particles from 
escaping, this can ensure the particles in upstream region have enough opportunities of scattering 
on double shocks to obtain more energy gains. In this way, we can obtain the maximum particle 
energy beyond lOMeV to investigate the energy spectral “break” between 1-lOMeV. 

In this simulation box, the continuous bulk flow from the middle of the box with opposite di¬ 
rections enter into the box along to the two boundaries of the simulation box. One bulk flow forms 
a CME shock in the left boundary, another bulk flow forms a bow shock on the right boundary. 
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With the shock formations, the particles are accelerated from each shock front. With the shock 
propagations forwarding to the middle of the box, the precursor regions are mixed. This interac¬ 
tion of the precursor regions ahead of two shocks make the energetic particles either re-accelerate 
or de-accelerate at the end of the simulation. It is this double effect of the energetic particles is 
responsible to the formation of the energy spectral “break” at the range of IMeV-lOMeV. 


In technically, the hybrid sim ulation method solves t he equation exp licitly for particle mo¬ 
tions in an electromagnetic plasma (IGiacalone et al.lll993l:lGuo et al.ll2013t) . But the Monte Carlo 
method applies the scattering law for the particles diffusive processes. According to the particles 
mean free path (mfp) equal to the local velocity times the scattering time (i.e. X = v/ • x), we deter¬ 
mine the scattering probabilities of particles based on the rate of the time step over the scattering 
time (i.e. r| = dt/x). So we can chose the numbers of the particles at a certain density in each grid 
to diffuse with their random pitch angle deflections. By means of these diffusive processes, the 
particles in upstream region transfer their kinetic energy into their random thermal energy in the 
downstream region. Then the minority of these random thermal particles can diffuse back shock 
front by multiple scattering cycles to obtain more energy gains to produce the superthermal energy 
“tail”. 


The simulation parameters mainly include upstream bulk speeds of m0i=- 0.6, and m 02=0.6, the 
total size of the box Xfnax=600, total simulation time rma;c=2400, the number of grids «x=1200, the 
initial density per grid no=360, the constant of the scattering time Xo=25/30, initial thermal velocity 
0)0 = 0.02, and time step dt=l/15. The above dimensionless values are all scaled by a group of the 
standard scaled factors: x^ca/e=2000/?g/600, t^(;a/e=63072400, and M^(;a/e=800kms“V0.6, where the 
Re is the Earth radius. The total number of the particles in simulation box amounts to more than 
1,000,000 particles. 


3. Results 
3.1. Shock Evolution 

The total simulation time rma;c=2400 are divided into 10 periods of time represented by Q=l, 
2, 3... 10 sequently. Fig J2| shows a group of bulk flow velocity snapshots in sequence of Q=3,5,7 
and density snapshots in sequence of Q=2, 4, 8. In the top panel, the bulk flow velocity with 
initial uOi and MO 2 in the upstream region evolute into the downstream with bulk flow speed of 
zero. From the Q=3, 5, to 7, the downstream region extend from the both boundaries to the center 
of the simulation box. Two shock fronts propagate from the both boundaries to the center of the 
box with opposite shock velocities v^hi and symmetrically. In front of the shocks, the velocity 

proflles have gradual slopes, which represent the shock precursors caused by the energetic particles 









Velocity Snapshot (Q=3) 


Velocity Snapshot (Q=5) 


Velocity Snapshot (Q=7) 



Grids Grids 




Fig. 2.— The top panel represents snapshots of the bulk flow velocity at the periods of simulation time 
Q=3, 5, 7. The low panel represents snapshots of the flow density at the periods of simulation time Q=4, 6, 

8 . 

back-reaction on the shock. When the two shock fronts approach more and more close, the two 
precursors overlap and have an interaction between the energetic particles extracted from the both 
downstream regions. This interaction process may vary the regular energy spectrum. In the low 
panel, the bulk flow density with an initial number dennsity no per grid evolute into downstream 
region with a higher density. The downstream region show a higher density with several times of 
that in the upstream region. Noted that both downstream regions have little lower densities in time 
period of Q=8 than Q=6 and 4, because our simulation model applying an adapter density with the 
simulation time. For the simply computation, we apply a temporally density with a reduction of 
density dn=16 in each period of time Q (i.e. adapter density reduce 1 particle from initial density 
no per 15 simulation time units). Similarly, the density profiles can also produce an interaction on 
both precursor regions when the two shocks become close enough. This may induce the precursor 
energy spectrum with the same spectrum in the downstream region at the end of the simulation. 

Fig j3] shows a group of the bulk flow velocity profiles at the periods of Q=3, 5, 7 and density 
profiles at the periods of Q=4, 6, 8. Top panel shows a series of the bulk flow velocity profiles 
in the position with the time. In the laboratory reference frame, the bulk flow velocity in both of 
the downstream regions are equal to zero, the bulk flow velocities in two upstream regions show 
an opposite movement in the middle of simulation box. The solid lines in both upstream regions 
represent the precursor position, respectively. The top right mesh plot shows the two precursor 
regions overlapping with time. This overlapping processes play two roles on the energetic particles: 
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Velocity profile (Q=3) Velocity profile (Q=5) Velocity profile (Q=7) 



Fig. 3.— The top panel represents profiles of the bulk flow velocity at the periods of simulation time 
Q=3,5,7. The low panel represents profiles of the bulk flow density at the periods of simulation time Q=4,6,8. 

(i) shock 2 with a positive precursor bulk velocity decelerate the energetic particles jumped from 
shock 1 with a negative precursor bulk velocity, (ii) In the contrary, shock 1 with a negative 
precursor bulk flow velocity re-accelerate the energetic particles dropped from shock 2 with a 
positive precursor bulk flow velocity. Lower panel shows a series of density profiles in the position 
with the time. The three mesh plots indicate the higher density in the downstream regions shorten 
the upstream regions with a lower density from the periods of Q=4, 6, to 8. The right plot shows 
two precursor regions are mixed together with a hybrid density. 


3.2. Particle Acceleration 

FigjU shows a group of particle acceleration processes in periods of Q=5, 6, 7, 8, 9, and 10. 
There are a part of particles extracted from the total simulation box in corresponding plot. In each 
pot, two triangle shadow areas represent the two shock fronts. Some of particles rotating in each 
shock front trace their trajectories with energy gains from the lower velocity to the higher velocity 
due to multiple scattering cycles on each shock front. Another particles in the downstream regions 
indicate they have no energy gains without velocity increasing. Maximum particle velocities Vmax 
denoted with values for 31.3898, 33.5277, 35.4340, 37.3926, 37.5091, and 36.4455 are calculated 
in periods of time Q=5, 6, 7, 8, 9, and 10, respectively. These maximum velocities show that these 
particles keep accelerating with the time to achieve an saturation at a certain time, then they exhibit 
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Particle Acceleration (Q=9) 



Fig. 4 .- This group of plots represent particles acceleration in double shocks at the periods of simulation time Q=5, 6, 7, 8, 9, and 10. These curves in each plot represent particle 

trajectories. These trajectories show local velocity evolutions of those particles with position and time. 


deceleration processes at the end of the simulation. In the periods of time between Q=8 and Q=9, 
two opposite shock fronts approach more close enough, it make the energetic particles produced 
by one shock cross into another shock each other. Simultaneously, these two precursor regions can 
make the energetic particles mixing in a hybrid precursor region. Because there exist two different 
modifications of the precursor bulk flow velocities, the energetic particles produced by one shock 
region may induce either deceleration or re-acceleration in another shock region. These mixture 
may break a smooth single energy spectral power law followed by a single shock model. 
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Fl^. 5 .- In top left plot, the solid line represents the proton spectrum with a “break” power law obtained from the downstream region at the end of simulation, the dashed line represents 

the thermal distribution in downstream region. In top right plot, the red line shows energy spectrum in the precursor region; the blue line shows energy spectrum in the downstream region. The 
middle left plot is the higher energy part of the top left plot. The pink line and cyan line show an energy spectral “break” at ~5MeV. The middle right plot also indicates double energy spectral 
indices. Cyan curve fitted the simulation data shows one index in lower energy range is less than 2, and another index in higher energy range is more than 2. The low left plot shows the energy 
spectra in downstream region purely with comparison of observed spectrum from the multiple spacecraft. The flux in left plot is statistic in per second, it is consist with the flux of the right plot in 
period of time 44 hours. 


3.3. Energy Spectra 

Fig j5] shows a group of proton spectra calculated from simulation box at the end of the sim¬ 
ulation. Each plot shows an energy spectral “break” at ~5MeV. The blue curve in top left plot 
represents the energy spectrum with a Maxiwellian peak at a few keV and with an extended power 
law energy spectrar‘tail” in downstream region. The energy spectral “break” occurs at the energy 
of ~ 5MeV. The dashed line indicates the potential Maxwellian distribution in the shocked down¬ 
stream. The top right plot represents two energy spectra in upstream and downstream regions, 
respectively. The pink color curve denotes the energy spectrum in upstream region and the blue 
color curve denotes the energy spectrum in downstream region. From the energy range of ~ IMeV 
to ~ lOMeV, the energy spectrum in upstream region with the same shape of the energy spectrum 
in downstream. Both of them appear energy spectral “breaks” at ~ 5MeV. For convenience to 
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study the energy spectral “break”, we extract a part of energy spectrum from the total downstream 
energy spectrum as shown in the middle left plot. The pink line represents the harder energy spec¬ 
trum with an index less than 2, and the cyan line represents the softer energy spectrum with an 
index more than 2. The “break” indicate that there exist a double power law spectrum. The middle 
right plot gives an energy spectral shape in the flux of ■F{E) representation, in which the lower 
energy spectrum exhibits a descent slope shape, and the higher energy spectrum exhibits a dropped 
slope shape. The fitting curve clearly shows that there exist a double power-law with a “break” 
at ~ 5MeV. The low panel presents a comparison of the simulated energy spectrum with the ob¬ 
served energy spectrum. The low left plot is the simulated energy spectrum in the downstream 
region purely; the lo w right plot is the observed energy spectrum in the downstream region too 
( Mewaldt et al.ll2008 1. The difference of the flux between two energy spectra is caused by different 
integrating time. The simulated flux is integrated in one second, but the observed flux is the sum 
of the period from 2006 Dec 13, 02:00 to Dec 14, 22:00, which is equivalent to 1.584x lO^s. So 
the simulated and the observed fluxes are well agreement in the normalized time. 


4. Summaries and Conclusions 

In summary, we do simulation of the converged two shocks for obtaining the proton spec¬ 
trum directly. Comparing with the observed energy spectrum from multiple spacecraft. Simulated 
energy spectrum exhibits the consist “break” at the energy ~5MeV. Although the numerical com¬ 
putation is very expensive, we obtain the highest energy spectral “tail” up to ~10MeV. We have 
updated the results of the maximum particle energy, which can not predict the energy spectral 
“break” in a single shock model. Why the previous efforts can not predict the energy spectral 
“break” ? There would be two reasons: (i) the single shock model can not provide enough en¬ 
ergy to transfer the energetic particles to the highest energy up to ~10MeV. (ii) The shortage of 
the interaction mechanism make it impossible to re-accelerate or decelerate the energetic particles. 
However, the converged shocks model can satisfy these two essential conditions. Firstly, the con¬ 
verged shocks model can provide enough energy injection than a single shock model to produce 
the highest energy spectral “tail”. Secondly, it can also provide an interaction mechanism to break 
the standard single power law energy spectrum formed by a single shock model. In converged 
two shocks, those energetic particles produced in each individual shock can mix together into the 
hybrid precursor region when two shocks approach more and more close. It is just these oppo¬ 
site precursor bulk flow velocities present the variation of the energetic particles distribution and 
produce an energy spectrar‘break” at the higher energy tail. 
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